#include "LinearInterpolation_PGR.h"
#include <thrust\transform.h>
#include <thrust\device_ptr.h>
#include <thrust\fill.h>


LinearInterpolation_PGR::LinearInterpolation_PGR(int width, int height){
	this->width = width;
	this->height = height;

	cudaMalloc(&device_in, width * height * sizeof(float));
	cudaMalloc(&device_out, width * height * sizeof(float));
	cpu_out = new float[width * height];
}


LinearInterpolation_PGR::~LinearInterpolation_PGR(){
	cudaFree(device_in);
	cudaFree(device_out);
	delete[] cpu_out;
}



__global__ void interpolateKernel(float* in, float* out, int width, int height, const int RANGE){
	int x = (blockIdx.x * blockDim.x + threadIdx.x);
	int y = (blockIdx.y * blockDim.y + threadIdx.y);

	
	int x_range = 0;
	int y_range = 0;
	float x_diff = 0.0;
	float y_diff = 0.0;
	float origin = in[x + y * width];

	if(origin > 50){
		out[x + y * width] = origin;
		
		//Find nearest valid data in x-direction
		for(int i = 1; x + i < width && i < RANGE; i++){
			if(in[(x + i) + y * width] > 50.0){
				x_range = i;
				x_diff = (in[(x + i) + y * width] - origin) / x_range;
				break;
			}
		}

		//valid data within x ~ +4
		if(x_range > 0){
			for(int i = 1; i < x_range; i++)
				out[(x + i) + y * width] = x_diff * i + origin;
		}


		//Find nearest valid data in y-direction
		for(int i = 1; y + i < height && i < RANGE; i++){
			if(in[x + (y + i) * width] > 50.0){
				y_range = i;
				y_diff = (in[x + (y + i) * width] - in[x + y * width]) / y_range;
				break;
			}
		}

		//valid data within y ~ +4
		if(y_range > 0){
			for(int i = 1; i < y_range; i++)
				out[x + (y + i) * width] = y_diff * i + origin;
		}

		////fill empty diagonal xy-direction
		//for(int i = 1; i < y_range; i++)
		//	for(int j = 1; j < x_range; j++)
		//		if(in[(x + j) + (y + i) * width] < 50.0)
		//			out[(x + j) + (y + i) * width] = 
		//				((x_diff * j + origin) + (y_diff * i + origin)) / 2.0;

		//fill empty diagonal xy-direction
		float step = (float)x_range/(float)y_range;
		float x_steprange = (float)x_range;
		for(int i = 1; i < y_range; i++)
			for(int j = 1; j < (int)(x_steprange+0.5); j++){
				if(in[(x + j) + (y + i) * width] < 50.0)
					out[(x + j) + (y + i) * width] = 
						((x_diff * j + origin) + (y_diff * i + origin)) / 2.0;
				x_steprange -= step;
			}
	}	
}

float* LinearInterpolation_PGR::interpolate(float3* in){
	thrust::transform(
		thrust::device_ptr<float3>(in),
		thrust::device_ptr<float3>(in) + width * height,
		thrust::device_ptr<float>(device_in),
		float3Tofloat());

	thrust::fill(
		thrust::device_ptr<float>(device_out),
		thrust::device_ptr<float>(device_out) + width * height,
		0.0);

	interpolateKernel<<<dim3(width / 32, height / 24), dim3(32, 24)>>>
		(device_in, device_out, width, height, 5);

	interpolateKernel<<<dim3(width / 32, height / 24), dim3(32, 24)>>>
		(device_out, device_in, width, height, 5);
	
	interpolateKernel<<<dim3(width / 32, height / 24), dim3(32, 24)>>>
		(device_in, device_out, width, height, 5);


	/*postProcessKernel<<<dim3(width / 32, height / 24), dim3(32, 24)>>>
		(device_out, width, height);*/

	return device_out;
}